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Locating apparent horizons is not only important for a complete understanding of numerically 
generated spacetimes, but it may also be a crucial component of the technique for evolving black- 
hole spacetimes accurately. A scheme proposed by Libson et at, based on expanding the location 
of the apparent horizon in terms of symmetric trace-free tensors, seems very promising for use with 
three-dimensional numerical data sets. In this paper, we generalize this scheme and perform a 
number of code tests to fully calibrate its behavior in black-hole spacetimes similar to those we 
expect to encounter in solving the binary black-hole coalescence problem. An important aspect of 
the generalization is that we can compute the symmetric trace-free tensor expansion to any order. 
This enables us to determine how far we must carry the expansion to achieve results of a desired 
accuracy. To accomplish this generalization, we describe a new and very convenient set of recurrence 
relations which apply to symmetric trace- free tensors. 
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I. INTRODUCTION 



A major goal of numerical relativity is to simulate the coalescence of an orbiting pair of black holes. In studying such 
systems, we will be interested in determining many quantities: the energy and momentum radiated, the associated 
waveforms, the total angular momentum of the system, etc. In addition to these common physical quantities, we will 
also want to understand the causal structure of the spacetime. Not only will this give us a more complete picture 
of the dynamics, but it also seems that tracking the causal structure may prove to be a crucial step in successfully 
evolving black- hole spacetimes [jj]. Knowing which events are inside black holes may allow them to be excised from 
the computational domain, thereby avoiding numerical difficulties that have plagued black hole evolutions. Ideally, 
we would like to be able to track all of the event horizons in a given spacetime. However, this is not possible during 
the evolution: event horizons can only be reconstructed after the evolution is complete. Instead, apparent horizons 
can be located on each individual space-like hypersurface during an evolution. Since apparent horizons must lie inside 
event horizons and asymptote towards them as the system settles down, they provide much of the desired causal 
information. They can also be used to define regions that can be excised from a computation. 

Various methods exist for locating apparent horizons. In practice, one searches for marginally outer-trapped 
surfaces (MOTS) J|. The apparent horizon is the outermost such surface. For axisymmetric problems, shooting 
methods have been the most widely used for locating apparent horizons, although decomposition into orthogonal 
polynomials [Q, the solution of elliptic boundary- value problems ||||, and the use of curvature flows jl0| have 
been used. Unfortunately, shooting methods do not generalize to three-dimensional spatial slices. The first general 
apparent-horizon finders were based on a spherical-harmonic decomposition of the MOTS [pT[-^3[. In this approach, 
each coefficient in the spherical-harmonic expansion is determined, iteratively, by performing an integral over a 
complicated function that characterizes the MOTS. The MOTS equation can also be posed as an elliptic equation 
for a function that parametrically specifies the location of the MOTS |l4|] . Curvature flow methods are also certainly 
applicable in the general case of a three-dimensional spatial hypersurface. 
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Recently, a variant of the spherical-harmonic decomposition method has been proposed by Libson et al. []15[ . This 
approach is conceptually appealing, having two particularly nice features. First, the coefficients in the spherical- 
harmonic expansion are determined by a minimization procedure, eliminating the need to perform surface integrals. 
Second, Libson et al. have proposed the use of symmetric trace-free (STF) tensors for parametrically representing the 
MOTS. This latter feature is particularly appealing when Cartesian coordinates are used on the three-dimensional 
hypersurface. 

In this paper, we will review the method proposed by Libson et al. and describe how we generalize the method 
by extending the expansion in STF tensors to arbitrary order. Because one is always using a truncated expansion, 
it is important to understand clearly the behavior of the apparent-horizon finder when the maximum order of the 
expansion is varied. As we will see, the number of points where the MOTS is determined will also affect the behavior 
of the apparent-horizon finder. We have examined both of these effects in detail. 

The paper is organized as follows: In § [II] we outline the method and basic equations, and in § III we explain our 
numerical implementation. In § IV we carefully discuss results from various test calculations, and in § M we briefly 
summarize the most important results. All technical details areprovided in the appendices. Appendix ET contains a 
number of useful equations relating to STF tensors. Appendix || describes the storage for arbitrary rank tensors. In 
Appendix ^ we derive recurrence relations for STF tensors. Finally, in Appendix ^ we derive an expression for the 
area element on the apparent horizon. 



II. METHOD AND BASIC EQUATIONS 

A MOTS is a closed two-surface embedded in a three-dimensional spatial hypersurface and, therefore, can be 
defined as a level surface of some scalar function r. If we use Cartesian coordinates on the hypersurface then we can 
parametrically define the level surface as 



r{x,y,z) = yJSijtx* - C*)(xi - Ci) - /(#,</>) = 0. (1) 

Here, x 1 are Cartesian coordinates, C* is a location inside the r = surface, and 9 and cf> are polar coordinates 
centered on C l . The function f(9, cj>) then measures the coordinate distance between C l and the r = surface in the 
direction (9, cj>). The outward pointing unit normal on the r = surface is 

Si = XdiT, (2) 

where A is the normalization factor 

A=[ s «(giT)(a,T)]~* (3) 

and gij is the metric on the spatial hypersurface. The expansion O of an outgoing null-bundle can now be written 

9 = Di S l + Kij S l S j -Kl (4) 

where Kij is the extrinsic curvature and Di is the covariant derivative operator associated with g^. Note that the 
term DiS 1 involves both first and second derivatives of r and hence of f(9, </>). Writing these out explicitly yields 



6 = - S l Si) - n inj ) - Xdidjf - SkT^j - K^J . (5) 

Here Sij is the Kronecker delta, n l is the unit vector in the (9, <fi) direction, and are the connection coefficients 
associated with . One must be careful in deriving this equation since there are effectively two metrics being used: 
the spatial metric g^ and the Kronecker delta Sij used in computing coordinate distances. In particular, we note that 

ryl f~<i 

(6) 



y/6 jk (xi - Ci)(x k - C k ) 

m = Si 3 n j (Sijrtrj = 1) (7) 

Si - X(m - dif) (8) 

S^g^Sj f/,,.S'.V 1; (9) 
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Our goal is now to find a function f(6, 4>) such that the r = surface is a MOTS, i.e. that the expansion (g) 
vanishes on that surface. In practice, instead of making the expansion vanish, we can evaluate it at a number Nq of 
points on the surface r = and look for a function / such that 



S(N e ) = L W a &i (10) 

a=l 

vanishes. If ([!(]) vanishes for arbitrary weights W a , then in the limit Nq — ► oo (so as to completely cover the surface) 
we are guaranteed of having located a MOTS. 

Our strategy will be to expand f(9, <fi) in terms of multipolc moments and to search for a minimum in S. The 
sum ( |l0| ) then depends on the corresponding expansion coefficients, which can be varied until the sum assumes a 
minimum. If this minimum comes arbitrarily close to zero an apparent horizon has been found. Thus, the problem 
has been reduced to a multidimensional minimization, for which standard methods can be used. 

An obvious choice of basis functions for the expansion of f{9, 4>) are the spherical harmonics, 

l i 

/(M)=E E F em Y lm (9,<t>), (11) 

1=0 m = -t 

where the expansion is truncated at order L. However, since we have to take up to second derivatives with respect to 
Cartesian coordinates, an expansion in terms of STF tensors 

L 

m<t>) = Y l J 7 K t NK t (12) 

1=0 

turns out to be a better choice. In the following we adopt the notation of where additional details of this 
formalism can be found. In particular, repeated indices will always be summed over. The subscript Kt denotes a 
multi- index of length £, and Nx e is the vector product of I unit vectors nf. 

N Ke = n kl n k2 ■ --nkf (13) 

In (|l2|), these are contracted with the STF tensors J-K e (of rank I). These are the location-independent expansion 
coefficients equivalent to the F lm in (0). Note that a STF tensor of rank I has 2£+ 1 independent components, just 
like the spherical harmonics. The relationship between ([ll]) and ( |l2|) can be seen even more clearly by choosing the 
basis for the STF tensors as defined in Ref. Jl6| (see Appendix^). In terms of these, !Fr 1 can be written 

i 

T Kl = J2 F lm y l £, (14) 

where the F lm are the same as in (pd[). The J^k™ also provide a relation between the spherical harmonics and the 

N Ke 

Y em (9,<p)=y e £N K M<f>)- (15) 

Inserting this into ([ll]) and using ( [Ti]) immediately yields (p"2|). 

Note that because the coefficients are location independent, derivatives of / can be calculated from derivatives 
of Nk c 

L 

d i f{e, ( ( > ) = J2^KA N K t (i6) 

£=0 

and similarly for second derivatives. 
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III. NUMERICAL IMPLEMENTATION 



Our code is designed in such a way that it can find a MOTS to arbitrary order L in the multipole expansion. On 
input we therefore have to specify the order L. Also, we have to specify the number of points Nq on the surface at 
which the expansion (||) (and, of course the sum for S) are evaluated. These points must be distributed somehow 
over the surface. Currently, they are distributed equally in <j) and cosf? on the unit sphere, but different choices could 
easily be made. Results for different values of L and Nq will be presented in the next section. 

Next the tensors Nx t , their first and second derivatives, as well as the basis STF tensors y^ have to be initialized. 
The latter are independent of location, so that we need to calculate them only once. The Nk { do, however, depend 
on the direction (#, </>), and therefore have to be calculated once for every point on the surface. In the code, we define 
arrays of length Nq to store Nx e and its derivatives. 

Since all of the STF tensors are completely symmetric, we can store the independent components in a very elegant 
and efficient way. This is explained in detail in Appendix [b]. In Appendix ^ we present recurrence relations that 
allow for a very efficient initialization of these objects. 

The MOTS search is started with a set of trial expansion coefficients F lm . These are then contracted with the 
basis STF tensors y^, which yields the Tk^ Since these quantities are independent of location, this needs to be 
done only once per iteration step. 

The are then contracted with Nx e and its derivatives to find /, dif and didjf for each direction (9, <fi). Once 
/ is known we can construct the coordinate location 

x l = fn l + C l (17) 

of the trial surface (r = 0) at each of the Nq points on the surface. For each point, we read in g^, and 1^. These 
can either be numerically evolved quantities or, for the test purposes in this paper, analytical values. Eq. (|J) now 
yields the expansion O for this location on the trial surface. Repeating these steps for every point Nq we can finally 
construct the sum S. Currently, we choose the weights W a in ( ^f]| ) based on the proper area element ( D13| ) defined 
in Appendix [d]. Thus, equation ([Io|) is an approximation to the mean square of the expansion: 

S = £<d 2 d 2 (j (18) 

Any multidimensional minimization routine can now be used to vary the F im until S has assumed a minimum. So 
far we have found best results with Powell's method |L7|], although it is likely that a method that uses derivatives 
with respect to the expansion coefficients F will be significantly faster, especially when the initial guess is already 
close to the final answer. We hope to explore this in the future. 

Once a minimum has been found we shift the center of the black hole C % according to the dipole moment (i.e. we 
choose C % so that the I = 1 moments vanish). We then repeat the MOTS search until the 1=1 moments stay below 
a predetermined maximum value. This procedure enables us to locate apparent horizons even when the initial guess 
is very poor (see next section) and should allow us to follow black holes that move through a numerical grid. 



IV. TESTS 



A. Schwarzschild 



An obvious test for the apparent-horizon finder is the Schwarzschild spacetime. Since the MOTS is spherically 
symmetric, it can be described with the monopole term alone. 

This test strongly demonstrated how well the shifting of the center C % according to the £ = 1 moments works. 
The code was able to locate the MOTS accurately even when the initial guess was completely disjoint from the true 
horizon. The code worked equally well when we located the black hole away from the origin of the coordinate system. 
In all cases the sum S, as well as all expansion coefficients F lm with I > 0, vanished to whatever tolerance we specified. 



B. Two black holes 



A spacetime containing two black holes has multiple MOTS, some of which can be highly distorted. Such a spacetime 
provides a much stronger test for the apparent-horizon finder. 

A metric for two time-symmetric black holes can be written in the conformally-flat form 
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ds 2 = V> 4 (dx 2 + dy 2 + dz 2 ), 



(19) 



where the conformal factor iji is given by 



2ri 2r 2 



and n and r 2 are 

ri = (x 2 + y 2 + (z + z j 2 ) 1/2 

r 2 = {x 2 +y 2 + (z~z ) 2 ) 1/2 . (21) 

Here M is the mass of the individual black holes, and zo is their coordinate distance from the origin of the coordinate 
system. 

Note that the singularities in ( |20| ) can be removed by adding matter sources (see Ref. |(| ) . Since this is advantageous 
in a numerical application and does not change the external metric we have implemented this form of the equations. 

The causal structure of this spacetime has been investigated in detail by Bishop Q . The MOTS can be found by 
using a shooting method to solve a set of coupled differential equations to high accuracy. This provides us with a 
solution that we can check the STF based apparent-horizon finder against. In the following we will refer to these 
solutions as the "true horizons" . 

In general there will be MOTS around the individual holes. If the holes are close enough, a pair of encompassing 
MOTS will also appear (see Ref. 0] for a careful discussion). According to Cadez these encompassing MOTS first 
appear at as separation of zq = 0.765. For separations close to the critical separation (i.e. zq < 0.765) they will be 
strongly distorted. 

Note that this situation is quite similar to what we expect in a binary black-hole evolution. Having that application 
in mind, it provides a strong test for our apparent-horizon finder and it can help us to decide to which order L we 
need to expand and at how many points Nq we need to evaluate the expansion O in order to accurately locate the 
an encompassing MOTS. 

In Fig. |l| we plot the estimated location of the MOTS based on expansions to order L — 2, 4, 6 and 8 for the case 
zq = 0.74. (By symmetry, only even L can contribute since C l — > 0). The "true horizon" is fairly distorted, causing 
the lower order expansions to perform very poorly. 

Note that all the lower order expansions find a location for the MOTS outside of the true horizon. This behavior 
could have severe consequences in numerical evolution codes that use "apparent horizon boundary conditions" (!]] and 
that ignore the causally disconnected region inside a MOTS. If we were to use one of the lower-order MOTS solutions 
for this purpose we could be ignoring a region that is not causally disconnected. 

This test clearly demonstrates that high order expansion is absolutely necessary for the detection of highly distorted 
horizons. On the other hand it also demonstrates that high order expansion is very expensive: increasing the order 
of the expansion by 2 increases the CPU time by roughly a factor of 3, ranging from several seconds for L — 2 to 
several minutes for L = 8 (on a serial computer). Details will depend on the particular numerical implementation as 
well as parameters associated with a giv en minimization routine. However, since we are searching for minima in an 



(L + l) 2 -dimensional space (see eq. (B9)) the required CPU time will always be a steep function of L. 

Another important factor for both the accuracy and the CPU time is the number of points iVe at which the 
expansion O is evaluated. In Fig. ^ we show results for zq = 0.6 using different numbers of points Nq = ng x n$, 
where ng is the number of points in the ^-direction and n^ in the 0-direction. 

Since the F lm up to order L have (L + l) 2 independent components, we will need at least (L + l) 2 points. From 
Fig. H it is obvious that 8x8 points are not enough for an expansion to order L = 8: the result is worse than for a 
lower order expansion. However, it can also be seen that increasing the number of points beyond this minimum can 
drastically increase the accuracy. 

As a next test we plot in Fig. || the integral (|L8|) as a function of zq for different expansion orders. For all these 
calculations we start with an initial guess close to where we expect a common MOTS. For zq < 0.765, when the two 
black holes have a common MOTS, we would therefore expect this sum to vanish, if we could resolve the MOTS 
arbitrarily well. For values of zq larger than 0.765 the minimization routine will find a nonzero minimum - however 
we expect these minima to go to zero as we approach zq — 0.765. 

Since we use an expansion to finite order we cannot resolve the MOTS arbitrarily well. This means that, while we 
will find a minimum, it will typically be different from zero even for zq < 0.765. The value will, again, depend on 
a number of parameters, but primarily on the order of the expansion L. This can be seen very clearly in Fig. |3|. In 
particular for the lower order expansions the expected drop at zq = 0.765 cannot be detected at all - a significant 
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FIG. 1. The estimated location of the MOTS found with expansions taken to order L — 2, 4, 6 and 8 for zo = 0.74 (using 
iVe = 25 X 11 points). The "true location" was found independently by solving a set of coupled differential equations (see 
Ref. [n]). The dot marks the coordinate center of one of the two black holes. 
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FIG. 2. Relative errors Ar/r as a function of cos# for expansion to order L = 2, 4, 6 and 8 and for different values 
iVe = no x (zo = 0.6). 
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decrease can only be seen for 1 = 8. This demonstrates again that an early detection of a common MOTS will only 
be possible with a high order expansion. 

Also, this suggests that it is hardly possible to decide on the basis of the value of the sum S if a MOTS has 
been found. As a better test, we suggest checking whether O is negative everywhere on a surface just inside the 
approximate MOTS. This surface "just inside" can be determined very easily by reducing the monopole term F 00 by 
a small fraction. 



V. SUMMARY 



We have developed an apparent-horizon finder based on a multipole expansion to arbitrary order L. The primary 
application we have in mind is the numerical evolution of a binary black-hole system. In order to check the performance 
of the MOTS finder in a spacetime of similar structure we have performed careful tests using initial data for two time- 
symmetric black holes. 

From these tests it is evident that a reliable search for highly distorted MOTS requires high order expansion. On 
the other hand, using a high order expansion is very expensive and it is questionable if this will be affordable during 
a dynamical evolution. 

However, in the evolution of a binary black-hole system for example, it is desirable to detect a common MOTS as 
early as possible since the region interior to this surface no longer needs to be evolved. It may, therefore, be worthwhile 
searching for this common MOTS using a high order expansion. 

As a compromise, it is possible to use a low order expansion for nearly spherical MOTS (as will be the case for 
the MOTS around individual black holes or the common MOTS in the later phases of an evolution) and a high order 
expansion for highly distorted MOTS (as in the early phase of the common MOTS). Unfortunately, "nearly spherical" 
is a coordinate dependent concept in this context, and it is not clear how well the coordinates will behave in a binary 
black-hole evolution code. 

An obvious optimization of the code is to change to a minimization scheme that uses derivatives. This is a fairly 
tedious, though straightforward task. 
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APPENDIX A: USEFUL EXPRESSIONS FOR STF TENSORS 



The purpose of this appendix is to provide several useful expressions for spherical harmonics and STF tensors that 
have been omitted in the main text. All these formulae are taken from reference | p6[ . 
The spherical harmonics can be written 

Y em (6, cf>) = C* n e im ^P £m (cos 6>) for < m < i 

[(l-m)/2] 

= C tm (e i 'l > sm0) m a em i{cos8) e - m - 2 J for < m < i (Al) 

3=0 

= (-l)">y*M* for -I < m < 0. 

Here * denotes complex conjugation, \{l — m)/2] means the largest integer less than or equal to [l — m)/2, and 

c^=(-ir( 2 \ +1 fr m l) 1/2 (A2) 



47T (^ + to)!_ 

/m] = (-I)'" (2*-2j)! 

2*j\{l-j)\{l-m-2j)V 1 1 
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Note that the Cartesian components of a radial unit vector n can be written 

n x + in y = e 1 ^ sin 9, n z = cos 9. (A4) 
In terms of these, the spherical harmonics are 

Y tm (6,<t>)=y%N Kt (0,<t>), (A5) 

where the are defined as 

l(t-m)/2] 

y£ = C *» £ ^ mj ^k 1 +^ kl )---^ {km +iSf k J 

= (-l) m (^r')* for-^<m<0. 

5%^ is the Kronecker delta and the brackets A( kl ... kl ) denote complete symmetrization. Obviously, the y^ are 
completely symmetric by definition, and they are also completely trace-free. Moreover, for every £ the 2£+l different 
y^ are linearly independent and span a basis for the STF tensors of rank I. Any STF tensor of rank I can therefore 
be written as 

e 

T Kt = J2 Fem y e £- (A7) 

m=-t 

For T Kl to be real the coefficients have to satisfy F e ~ m = (-l)™(F^ m )*. 

APPENDIX B: STORAGE 

Since we want to allow for expansion to arbitrary order L, we have to store objects of the type N K( with multi- 
indices up to arbitrary length L. This can be accomplished by noting that all these objects are completely symmetric, 
i.e. 

N kl ... ki ... k] ... kl — N kl ... kj ... ki ... ki (Bl) 

for any pair of indices fc, and kj. In three dimensions these indices can only take the value x, y or z. It is therefore 
sufficient to specify the number of indices with value x, y and z to determine any element in Nx e uniquely. In the 
following we adopt the notation 

N( Nx )( Ny )( Nz) , (B2) 
where the N x , N y and N z are the total numbers of indices x, y and z. Obviously, the rank of the tensor is 

£ = N X + N V + N Z . (B3) 

A symmetric tensor of rank £ has 

+ = \(e + W + 2) (B4) 
independent elements, the number of independent components for all symmetric tensors of rank up to L is 

L 

2 ' " () 



J2 \{i + m + 2) = \{L+ 1)(L + 2)(L + 3). (B5) 



We can therefore store all symmetric tensors up to order L in an array of this length. Each element is uniquely 
determined by a combination of N x , N y and N z . 
Note that an element ^(N !C )(Ny)(N z ) appears 
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w (B6) 

times in a symmetric tensor. This weighting factor has to be taken into account when carrying out sums as in (p^). 

In our code we store Njc t , its first and second derivative and the STF tensor Tki in arrays of this kind. In addition 
to being symmetric the latter is completely trace-free, so that we could use even less storage. We decided not to do 
so, since this would complicate the code for little benefit. 

The next objects that need to be stored are the basis STF tensors y^l- For each I this requires storage of 21 + 1 
values. We therefore need an array of length 

L 1 1 

Y, ^ + m + W + 2) = ^(3L + 2)(L + 1)(L + 2)(L + 3). (B7) 

1=0 

Each element is uniquely labeled by a combination of N x , N y and N z together with an index m. 

As a further complication these tensors are complex. However, since = (—l) m (yj<^)* (see appendix A, 

equation ( |A6| )), we need to store only the non- negative m. This can be accomplished by storing the real parts in the 
storage for m > and the imaginary parts in the storage for m < 0. 

A sum as in ( |l4| ) can be written 

e 

J~K t = 2^ b *Kt 

m=-t 

i 

= F m y™ + J2 [F lm y l K t + (F em yp;y] (bs) 

m— 1 

e 

= F eo yg e + 2^2 [MF em MyiQ - %F em zy e £] 

m— 1 

(where we have assumed F l ~ m = (_i) m (^™)* so that Tk,_ is reai ) see appendix A). Storing the elements as described 
above, we can ag ain sum from —i to £, but we have to take into account a weight of 2 for m > and —2 for m < 0, 



as given by eqn. (BJ 



The last objects that have to be stored are the expansion coefficients F im . Since for each £ there are 2£ + l different 
elements, we have to store a total of 

L 

]T(2^+1) = (L + 1) 2 (B9) 

elements. Again, these elements will be complex and satisfy F l ~ m = so that we can again store the 

real parts in the storage for m > and the imaginary parts in the storage for m < 0. 



APPENDIX C: INITIALIZATION AND RECURRENCE RELATIONS 



Before a MOTS can be found Njc t , its derivatives and y^ 1 need to be ini tiali zed for various directions 



Although direct expressions could be used to do so (as for example equation (A6) in Appendix A), this would be 
extremely tedious and inefficient. Moreover, sums like ( |A6| ) are prone to cancellation error from adding terms of 
opposite sign. We have therefore derived recursive expressions that make the initialization much easier. 

Since the yf^ 1 are closely related to the spherical harmonics, we can start with the standard recurrence relation 



Y*" = yj f±L [VW-IY^- cos \j ie ~ 2 Z7 n2 y (i - 2)m ) ■ (CD 

Here we can use cos 9 = n z , insert v? x + r? y + n 2 z — 1 , and replace the spherical harmonics with equation ( |l5| ) , which 
yields 



= V W^? V2F^iyt^N Ke _ inz - y {£ y^ m N Kt _ 2 (nl + nl + nl)). (C2) 
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a 



n In m \ J/ ( N -)( N y)( N ^ N( - N ^ N y^ N ^ 

x • y • z • 



Adopting our notation (B2), we can rewrite this equation as 
(£-1)1 



21+1 

P — m 2 



(C3) 



N x \N y \(N z - 1)! V LJ/ (N x )(N y )(N,-i) 



(£-2)\ 



(e-2)r 



2£-3 
(£-2)\ 



[£-2)\ 



y 



(l-2)r. 
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Note that, because the are direction- independent, only specific terms in the STF tensor sums are allowed on 

the right-hand side of eq. ( |C2[ ). Only those terms that, together with the additional unit vectors appearing on the 
right-hand size of eq. ( |C2| ) , result in the same number of N x , N y , and N z that appear on its left-hand size are allowed. 
Finally, the direction- independence of the 3^™, means that we can drop the the -/V(Af x )(jv B )(JV,) terms (which contain 
all of the directional dependence) to get 
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This expression can be used to initialize the y^ with — I < m < £. For m = £ we can use eq. ( |A6| ), which for this 
case reduces to 
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Next we have to initialize the tensors Nx e (eq. (|l^)) as well as their first and second derivatives. Note that the 
partial derivative of a unit vector is 
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The partial derivative of Nk £ is therefore 

diN Kl = 9i(n kl nk 2 ■ --n^) = £n i ( kl N Kt _ 1 y 

Taking a second derivative yields, after some algebra, 
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didjN Ke = — {nijN Ki + ^n^^e-i) + rijn^ kt N k^)) + £(t - ^)n i i ke _ 1 N Ke _ 2 n ke)j . 
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Note that we can now construct Njct, diNK e and didjNK t from the ten different, totally-symmetric objects Nk ( , 
n i(kt^K e -i) an< i n i{ke- 1 NK e ^2 n ki i )j- Denoting any one of these with ^Aorjv HAy (in our notation (B2)) we find that 
they all satisfy the recurrence relation 
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N x + N y + N z 

The starting values for the different objects are 

Nx e n a ^ ke N Kei ) n a ( kl _ 1 N Kl _ 2 n kll ) b 

I = 1 

1 = 1 rii n al 

1 = 2 mrij \ijiaiTij + n a jrii) \(n ai n ho + n aj n bi ) ■ 

Once the se objects have been calculated up to order L, we can then construct the derivatives of Nk { using eqs. (C7) 
and (pi). 
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APPENDIX D: THE AREA ELEMENT ON THE APPARENT HORIZON 

Beginning with the line element for the three-dimensional spatial hypersurface 

ds 2 = g i:j dx l dx j , (Dl) 
we transform to polar coordinates centered around some point using 

The line element for the surface parameterized by r = f(6, </>) is obtained by substituting 

dr = % d0 + lk d * ^ 



into eqn. (D2). We define 



and 



which can be expanded to yield 
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O x = -7==(1 + n x 8 x f + n y d y f) - n x y/T^ z dJ, (D6) 

&y = Zhg» (1 + ^aj + _ ^^T^Jc^/, ( D7 ) 

v 1 

9 

e 2 = j=JL= {n x d x f + n y d y f) - y/T=^(l + n z d z f), (D8) 

$ z = ^ — == [n x {n x d y f - n y d x f) - n y ], (D9) 
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VI 3 »_ 

$* = J— [ % (n a ^/ - + n*], (D10) 



^ = J— (n x d y f - n y d x f). (Dll) 

v 1 

Finally, we can now rewrite the line element ( |Dl| ) as 

ds 2 = / 2 (e^dfl 2 + 26*$, sin 6W0# + sin 2 6>d0 2 ) , (D12) 
where we have lowered the indices on and <!> with the metric gij . From this we find the two-surface area element 



d 2 u = J(e i G i )(&$j) - (O k $ k ) 2 f 2 sindd6d(t>. (D13) 
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